#!pip install numpy==1.23.3 --user
# Import Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import plotly
import plotly.offline as pyoff
import plotly.graph_objs as go
import plotly.express as px
import chart_studio
import chart_studio.plotly as py
import calmap
import datetime
import tensorflow as tf
import os
import random
import re
import plotly.offline as pyoff
import plotly.graph_objs as go
import swifter
from scipy import signal
from datetime import date
from plotly.subplots import make_subplots
from itertools import cycle, product
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from scipy.stats import boxcox
from dateutil.parser import parse
from pmdarima.arima import auto_arima
from pmdarima.utils import diff_inv
from statsmodels.tsa.stattools import adfuller
from sklearn.model_selection import TimeSeriesSplit
from tensorflow.keras.layers import LSTM, Dense, BatchNormalization
from tensorflow.keras import Sequential
from tensorflow.keras.backend import clear_session
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator
from tensorflow.keras.initializers import *
from tensorflow.keras import optimizers
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, mean_absolute_error, r2_score
from sklearn.linear_model import LinearRegression
from scipy.special import boxcox1p, inv_boxcox1p
import matplotlib.patches as mpatches
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.model_selection import GridSearchCV
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from tensorflow.keras.optimizers import Adam
from IPython.display import HTML, display
from swifter import set_defaults
from xgboost import XGBRegressor
from xgboost import plot_importance, plot_tree
# Versões dos pacotes usados neste jupyter notebook
%reload_ext watermark
%watermark -a "Herikc Brecher" --iversions
Author: Herikc Brecher chart_studio: 1.1.0 tensorflow : 2.10.0 seaborn : 0.11.2 scipy : 1.5.4 calmap : 0.0.9 matplotlib : 3.5.1 re : 2.2.1 keras : 2.10.0 numpy : 1.20.0 swifter : 1.3.4 plotly : 5.6.0 pandas : 1.4.2
# Variaveis globais
SEED = 84796315
FEATURES = 7
EPOCHS = 1000
BATCH_SIZE = 1000
EXECUTE_GRID_SEARCH = False
# Configurando seeds
os.environ['PYTHONHASHSEED'] = str(SEED)
tf.random.set_seed(SEED)
np.random.seed(SEED)
random.seed(SEED)
# Exibindo toda tela
display(HTML('<style>.container { width:100% !important; }</style>'))
pd.options.plotting.backend = 'matplotlib'
# Configurando swifter
set_defaults(
npartitions = None,
dask_threshold = 1,
scheduler = "processes",
progress_bar = True,
progress_bar_desc = None,
allow_dask_on_strings = True,
force_parallel = True,
)
# Import dataset
dtOrders = pd.read_csv('../data/olist_orders_dataset.csv', encoding = 'utf8')
# Colunas do tipo data
dateColumns = ['order_purchase_timestamp', 'order_approved_at', 'order_delivered_carrier_date',\
'order_delivered_customer_date', 'order_estimated_delivery_date']
# Dataset de analise temporal
dtOrdersAdjusted = dtOrders.copy()
# Convertendo columas de data para date
for col in dateColumns:
dtOrdersAdjusted[col] = pd.to_datetime(dtOrdersAdjusted[col], format = '%Y-%m-%d %H:%M:%S')
# Dropando valores NA
dtOrdersAdjusted = dtOrdersAdjusted.dropna()
dtOrdersAdjusted.dtypes
order_id object customer_id object order_status object order_purchase_timestamp datetime64[ns] order_approved_at datetime64[ns] order_delivered_carrier_date datetime64[ns] order_delivered_customer_date datetime64[ns] order_estimated_delivery_date datetime64[ns] dtype: object
dtHistory = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date
start = dtHistory.min()
end = dtHistory.max()
idx = pd.date_range(start, end, normalize = True)
seriesOriginal = dtHistory.value_counts(sort = False).sort_index().reindex(idx, fill_value = 0)
dtHistory = pd.DataFrame(seriesOriginal).reset_index()
Principais outliers identificados:
dtHistory.rename(columns = {'index': 'Data', 'order_purchase_timestamp': 'Vendas'}, inplace = True)
dtHistory
| Data | Vendas | |
|---|---|---|
| 0 | 2016-09-15 | 1 |
| 1 | 2016-09-16 | 0 |
| 2 | 2016-09-17 | 0 |
| 3 | 2016-09-18 | 0 |
| 4 | 2016-09-19 | 0 |
| ... | ... | ... |
| 709 | 2018-08-25 | 69 |
| 710 | 2018-08-26 | 73 |
| 711 | 2018-08-27 | 66 |
| 712 | 2018-08-28 | 39 |
| 713 | 2018-08-29 | 11 |
714 rows × 2 columns
# Plot
# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
y = dtHistory['Vendas'])]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
yaxis = {'title': 'Vendas'},
title = 'Vendas por dia')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0, 800])
pyoff.iplot(fig)
# Remove outliers
seriesOriginal = seriesOriginal[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
pred_range = pd.date_range(datetime.date(2018, 8, 17), datetime.date(2018, 10, 17))
dtHistory = pd.DataFrame(seriesOriginal).reset_index()
dtHistory.rename(columns = {'index': 'Data', 'order_purchase_timestamp': 'Vendas'}, inplace = True)
# Plot
# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
y = dtHistory['Vendas'])]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
yaxis = {'title': 'Vendas'},
title = 'Vendas por Dia')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0, 800])
pyoff.iplot(fig)
# Plot
# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['Data'],
y = dtHistory['Vendas'],
showlegend = False),
go.Scatter(x = dtHistory['Data'],
y = -dtHistory['Vendas'],
fill='tonexty',
showlegend = False)]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
yaxis = {'title': 'Vendas'},
title = 'Vendas por Dia Expandida')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_traces(line_color = 'green')
fig.update_layout(yaxis_range = [-800 , 800])
pyoff.iplot(fig)
#Plot histórico de vendas por dia, mês e ano
fig, caxs = calmap.calendarplot(seriesOriginal, daylabels = 'MTWTFSS', fillcolor = 'grey',cmap = 'YlGn', fig_kws = dict(figsize = (18, 9)))
fig.suptitle('Histórico de Vendas', fontsize = 22)
fig.subplots_adjust(right = 0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.67])
fig.colorbar(caxs[0].get_children()[1], cax = cbar_ax)
plt.show()
# Criar grafico na estrutura STL 4 layers
def add_stl_plot(fig, res, legend):
axs = fig.get_axes()
# Nome de cada um dos subplots
comps = ['trend', 'seasonal', 'resid']
for ax, comp in zip(axs[1:], comps):
series = getattr(res, comp)
if comp == 'resid':
ax.plot(series, marker = 'o', linestyle = 'none')
else:
ax.plot(series)
ax.legend(legend, frameon = False)
# Gerar STL
stl = STL(seriesOriginal)
stl_res = stl.fit()
fig = stl_res.plot()
fig.set_size_inches((20, 12))
plt.show()
# Gerar STL não robusto e concatenar ao robusto
stl = STL(seriesOriginal, robust = True)
res_robust = stl.fit()
fig = res_robust.plot()
fig.set_size_inches((20, 12))
res_non_robust = STL(seriesOriginal, robust = False).fit()
add_stl_plot(fig, res_non_robust, ['Robusto', 'Não Robusto'])
# Additive Decomposition
additive_decomposition = seasonal_decompose(seriesOriginal, model = 'additive', period = FEATURES)
# Plot
additive_decomposition.plot()
plt.tight_layout(rect = [0, 0.03, 1, 0.95])
plt.rcParams["figure.figsize"] = (20, 10)
plt.show()
# Removendo tendência
detrended = signal.detrend(seriesOriginal)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
y = detrended)]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Vendas Diárias Removendo a Tendência')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
# Gerando STL para separar cada um dos componentes
stl = STL(seriesOriginal)
res = stl.fit()
# Separando seriesDeseasonal
seriesDeseasonal = res.observed - res.seasonal
# Separando boxcox
seriesBoxCox, lmbda = boxcox(seriesOriginal + 1)
seriesBoxCox = pd.Series(seriesBoxCox, index = seriesOriginal.index)
# Separando stationary
seriesResidual = seriesOriginal.diff(7).dropna()
Os testes abaixo concluiram:
O teste aceita a hipótese nula em que a série não é estácionária para os dados originais e deseasonal. Já para os dados residuais esses aceitaram a hipótese alternativa que os dados são estacionários.
ADF teste:
print("Os dados são estacionários?\n")
testResult = adfuller(seriesOriginal, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")
for p, v in testResult[4].items():
print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))
Os dados são estacionários? Valor Teste = -2.616 Valor de P: = 0.090 Valores Críticos: 1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança 5%: -2.866544718556839 - O dataset não é estacionário com 95% de confiança 10%: -2.5694353738653684 - O dataset é estacionário com 90% de confiança
print("Os dados deseasonal são estacionários?")
testResult = adfuller(seriesDeseasonal, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")
for p, v in testResult[4].items():
print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))
Os dados deseasonal são estacionários? Valor Teste = -2.536 Valor de P: = 0.107 Valores Críticos: 1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança 5%: -2.866544718556839 - O dataset não é estacionário com 95% de confiança 10%: -2.5694353738653684 - O dataset não é estacionário com 90% de confiança
print("Os dados boxcox são estacionários?")
testResult = adfuller(seriesBoxCox, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")
for p, v in testResult[4].items():
print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))
Os dados boxcox são estacionários? Valor Teste = -2.874 Valor de P: = 0.048 Valores Críticos: 1%: -3.441694608475642 - O dataset não é estacionário com 99% de confiança 5%: -2.866544718556839 - O dataset é estacionário com 95% de confiança 10%: -2.5694353738653684 - O dataset é estacionário com 90% de confiança
print("Os dados residuais são estacionários?")
testResult = adfuller(seriesResidual, autolag = 'AIC')
print("Valor Teste = {:.3f}".format(testResult[0]))
print("Valor de P: = {:.3f}".format(testResult[1]))
print("\nValores Críticos:")
for p, v in testResult[4].items():
print("\t{}: {} - O dataset {} é estacionário com {}% de confiança".format(p, v, "não" if v < testResult[0] else "", 100 - int(p[:-1])))
Os dados residuais são estacionários? Valor Teste = -6.802 Valor de P: = 0.000 Valores Críticos: 1%: -3.441834071558759 - O dataset é estacionário com 99% de confiança 5%: -2.8666061267054626 - O dataset é estacionário com 95% de confiança 10%: -2.569468095872659 - O dataset é estacionário com 90% de confiança
Toda a etapa de modelagem será considerada com 5 passos a frente de previsão.
# Controle de resultados de toda fase de modelagem
result = pd.DataFrame(columns = ['Algorithm', 'MSE', 'RMSE', 'MAE', 'MAPE', 'Mean_True_Value', 'Mean_Predict_Value'])
split_range = TimeSeriesSplit(n_splits = 8, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])
# Adiciona o registro ao dataset
def record(result, algorithm, mse = -1, rmse = -1, mae = -1, mape = -1, mrv = -1, mpv = -1, show = True):
new = pd.DataFrame(dict(Algorithm = algorithm, MSE = mse, RMSE = rmse, MAE = mae, MAPE = mape, Mean_True_Value = mrv,\
Mean_Predict_Value = mpv), index = [0])
result = pd.concat([result, new], ignore_index = True)
if show:
display(result)
return result
# Plot no formato de 4 layers, seguindo o STL para cada um dos modelos
def plot(index, pred, mse, title, fig = None, ax = None, ylabel = ''):
global seriesOriginal
empty_fig = fig is None
if empty_fig:
fig, ax = plt.subplots(figsize = (13, 6))
else:
ax.set_ylabel(ylabel)
ax.set_title(title)
patch_ = mpatches.Patch(color = 'white', label = f'MSE: {np.mean(mse):.1e}')
L1 = ax.legend(handles = [patch_], loc = 'upper left', fancybox = True, framealpha = 0.7, handlelength = 0)
ax.add_artist(L1)
sns.lineplot(x = seriesOriginal.index, y = seriesOriginal, label = 'Real', ax = ax)
sns.lineplot(x = index, y = pred, label = 'Previsto', ax = ax)
ax.axvline(x = index[0], color = 'red')
ax.legend(loc = 'upper right')
ax.set_ylim([0, 800])
if empty_fig:
plt.show()
else:
return fig
# Calculo para previsão e teste quando utilizado a série Original
def calcPredTestOriginal(train, pred, test):
return pred, test, 0
# Calculo para previsão e teste quando utilizado a série seriesDeseasonal
def calcPredTestseriesDeseasonal(train, pred, test):
# Removendo a sazonalidade da série e convertendo para o shape correto
last_seasonal = res.seasonal.reindex_like(train).tail(stl.period)
pred = pred + np.fromiter(cycle(last_seasonal), count = pred.shape[0], dtype = float)
test = test + res.seasonal.reindex_like(test)
return pred, test, 1
# Calculo para previsão e teste quando utilizado a série BoxCox
def calcPredTestBoxCox(train, pred, test):
# Reverdendo a normalização do boxcox
pred = inv_boxcox1p(pred, lmbda)
test = inv_boxcox1p(test, lmbda)
return pred, test, 2
# Calculo para previsão e teste quando utilizado a série Stationary
def calcPredTestStationary(train, pred, test):
# Calculando a diferença da sazonalidade
xi = seriesOriginal.reindex_like(train).tail(FEATURES)
totalLen = len(pred) + len(xi)
ix = pd.date_range(xi.index[0], periods = totalLen)
inv = diff_inv(pred, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)
inv = pd.Series(inv, index = ix, name = 'Vendas')
pred = inv.iloc[FEATURES:]
totalLen = len(test) + len(xi)
ix = pd.date_range(xi.index[0], periods = totalLen)
inv = diff_inv(test, FEATURES, xi = xi) + np.fromiter(cycle(xi), count = totalLen, dtype = float)
inv = pd.Series(inv, index = ix, name = 'Vendas')
test = inv.iloc[FEATURES:]
return pred, test, 3
def calculate_metrics_off(test, pred):
mse = mean_squared_error(test, pred, squared = True)
rmse = mean_squared_error(test, pred, squared = False)
mae = mean_absolute_error(test, pred)
mape = mean_absolute_percentage_error(test, pred) * 100
mtv = np.mean(test)
mpv = np.mean(pred)
return mse, rmse, mae, mape, mtv, mpv
def calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv):
mse_, rmse_, mae_, mape_, mtv_, mpv_ = calculate_metrics_off(test, pred)
mse.append(mse_)
rmse.append(rmse_)
mae.append(mae_)
mape.append(mape_)
mtv.append(mtv_)
mpv.append(mpv_)
return mse, rmse, mae, mape, mtv, mpv
# Report para Time Series Regressor, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportTSR(data, modelName, calcFunction):
global result
global figs
mse = []
rmse = []
mae = []
mape = []
mtv = []
mpv = []
title = modelName + ' - Time Series Regression'
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
gen = TimeseriesGenerator(train, train, FEATURES, batch_size = BATCH_SIZE)
X_train = gen[0][0]
y_train = gen[0][1]
lr = LinearRegression()
lr.fit(X_train, y_train)
X_pred = y_train[-FEATURES:].reshape(1,-1)
pred = np.empty(test.shape[0])
for i in range(len(pred)):
forecast = lr.predict(X_pred)
X_pred = np.delete(X_pred, 0, 1)
X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)
pred[i] = forecast
pred, test, indexPlot = calcFunction(train, pred, test)
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.tight_layout()
plt.close()
reportTSR(seriesOriginal.copy(), 'Original', calcPredTestOriginal)
reportTSR(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal)
reportTSR(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox)
reportTSR(seriesResidual.copy(), 'Stationary', calcPredTestStationary)
result
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | |
|---|---|---|---|---|---|---|---|
| 0 | Original - Time Series Regression | 5296.329792 | 60.784585 | 42.960631 | 25.553663 | 179.560484 | 168.008596 |
| 1 | Deseasonal - Time Series Regression | 5081.77776 | 58.208746 | 40.122247 | 23.488563 | 179.560484 | 164.670085 |
| 2 | BoxCox - Time Series Regression | 5328.602465 | 61.060147 | 43.254815 | 25.416806 | 179.560484 | 165.399796 |
| 3 | Stationary - Time Series Regression | 6057.828087 | 64.323443 | 46.538896 | 27.899678 | 179.560484 | 184.267746 |
# Função utilizada para o hypertuning de alpha, beta e gamma do Exponential Smoothing
def GSES(data, modelName, alpha, beta, gamma, damping_trend, calcFunction):
mse = []
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
try:
with catch_warnings():
filterwarnings('ignore')
ES = (
ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
.fit(smoothing_level = alpha, smoothing_trend = beta, smoothing_seasonal = gamma, method = 'ls', damping_trend = damping_trend)
)
pred = ES.forecast(test.shape[0])
pred, test, _ = calcFunction(train, pred, test)
mse.append(mean_squared_error(pred, test, squared = True))
except:
mse.append(-1)
return np.mean(mse)
# Função utilizada para o hypertuning de demais parâmetros do Exponential Smoothing
def GSESOPT(data, modelName, trend, season, periods, bias, method, calcFunction):
mse = []
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
try:
with catch_warnings():
filterwarnings('ignore')
ES = (
ExponentialSmoothing(train, trend = trend, seasonal = season, seasonal_periods = periods)
.fit(remove_bias = bias, method = method, optimized = True)
)
pred = ES.forecast(test.shape[0])
pred, test, _ = calcFunction(train, pred, test)
mse.append(mean_squared_error(pred, test, squared = True))
except:
mse.append(-1)
return np.mean(mse)
# Report para Exponential Smoothing, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportES(data, modelName, model_kwargs, fit_kwargs, calcFunction):
global result
global figs
mse = []
rmse = []
mae = []
mape = []
mtv = []
mpv = []
title = modelName + ' - Exponential Smoothing'
indexPlot = 0
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
ES = (
ExponentialSmoothing(train, trend = model_kwargs['trend'], seasonal = model_kwargs['seasonal'], seasonal_periods = FEATURES, damped_trend = model_kwargs['damped_trend'])
.fit(smoothing_level = fit_kwargs['smoothing_level'], smoothing_trend = fit_kwargs['smoothing_trend'],\
smoothing_seasonal = fit_kwargs['smoothing_seasonal'], method = fit_kwargs['method'], damping_trend = fit_kwargs['damping_trend'])
)
pred = ES.forecast(test.shape[0])
pred, test, indexPlot = calcFunction(train, pred, test)
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
# Função para gerar tabela de hypertuning ampla
def exp_smoothing_configs(seasonal = [None]):
models = list()
# Lista de argumentos
t_params = ['add', 'mul']
s_params = ['add', 'mul']
p_params = seasonal
r_params = [True, False]
method_params = ['L-BFGS-B' , 'TNC', 'SLSQP', 'Powell', 'trust-constr', 'bh', 'ls']
# Gerando lista de argumentos
for t in t_params:
for s in s_params:
for p in p_params:
for r in r_params:
for m in method_params:
cfg = [t, s, p, r, m]
models.append(cfg)
return models
%%time
# Treinamento do modelo
if EXECUTE_GRID_SEARCH:
# Gerando tabela de hypertunning
alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesOriginal.copy(), 'Original',\
x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestOriginal), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
%%time
# Se True irá treinar com a nova lista mais ampla (pode demorar)
if EXECUTE_GRID_SEARCH:
# Criando lista de argumentos ampla
params_ = exp_smoothing_configs([FEATURES])
# Convertendo dicionário de argumentos para dataframe
hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
# Executando hypertuning
hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesOriginal.copy(), 'Original',\
x.trend, x.season, x.periods, x.bias, x.method, calcPredTestOriginal),\
axis = 1)
# Busca em query pelos melhores argumentos
display(hyperparam_.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13, 6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.8, smoothing_seasonal = 0, method = 'ls', damping_trend = 0.8)
reportES(seriesOriginal.copy(), 'Original', model_kwargs, fit_kwargs, calcPredTestOriginal)
%%time
if EXECUTE_GRID_SEARCH:
alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesDeseasonal.copy(), 'seriesDeseasonal',\
x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestseriesDeseasonal), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
%%time
if EXECUTE_GRID_SEARCH:
params_ = exp_smoothing_configs([FEATURES])
hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesDeseasonal.copy(), 'seriesDeseasonal',\
x.trend, x.season, x.periods, x.bias, x.method, calcPredTestseriesDeseasonal),\
axis = 1)
display(hyperparam_.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.2, smoothing_seasonal = 0.5, method = 'ls', damping_trend = 0.8)
reportES(seriesDeseasonal.copy(), 'Deseasonal', model_kwargs, fit_kwargs, calcPredTestseriesDeseasonal)
%%time
if EXECUTE_GRID_SEARCH:
alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesBoxCox.copy(), 'BoxCox',\
x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestBoxCox), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
%%time
if EXECUTE_GRID_SEARCH:
params_ = exp_smoothing_configs([FEATURES])
hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesBoxCox.copy(), 'BoxCox',\
x.trend, x.season, x.periods, x.bias, x.method, calcPredTestBoxCox),\
axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.1, smoothing_trend = 0.7, smoothing_seasonal = 0.0, method = 'ls', damping_trend = 0.8)
reportES(seriesBoxCox.copy(), 'BoxCox', model_kwargs, fit_kwargs, calcPredTestBoxCox)
%%time
if EXECUTE_GRID_SEARCH:
alphas = betas = gammas = damping_trend = np.arange(1, step = 0.1)
hyperparam = pd.DataFrame(product(alphas, betas, gammas, damping_trend), columns = ['alpha', 'beta', 'gamma', 'damping_trend'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSES(seriesResidual.copy(), 'Stationary',\
x.alpha, x.beta, x.gamma, x.damping_trend, calcPredTestStationary), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
%%time
if EXECUTE_GRID_SEARCH:
params_ = exp_smoothing_configs([FEATURES])
hyperparam_ = pd.DataFrame(params_, columns = ['trend', 'season', 'periods', 'bias', 'method'])
hyperparam_['mse'] = hyperparam_.swifter.apply(lambda x: GSESOPT(seriesResidual.copy(), 'Stationary',\
x.trend, x.season, x.periods, x.bias, x.method, calcPredTestStationary),\
axis = 1)
display(hyperparam_.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
model_kwargs = dict(trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
fit_kwargs = dict(smoothing_level = 0.0, smoothing_trend = 0.2, smoothing_seasonal = 0.1, method = 'ls', damping_trend = 0.2)
reportES(seriesResidual.copy(), 'Stationary', model_kwargs, fit_kwargs, calcPredTestStationary)
result
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | |
|---|---|---|---|---|---|---|---|
| 0 | Original - Time Series Regression | 5296.329792 | 60.784585 | 42.960631 | 25.553663 | 179.560484 | 168.008596 |
| 1 | Deseasonal - Time Series Regression | 5081.77776 | 58.208746 | 40.122247 | 23.488563 | 179.560484 | 164.670085 |
| 2 | BoxCox - Time Series Regression | 5328.602465 | 61.060147 | 43.254815 | 25.416806 | 179.560484 | 165.399796 |
| 3 | Stationary - Time Series Regression | 6057.828087 | 64.323443 | 46.538896 | 27.899678 | 179.560484 | 184.267746 |
| 4 | Original - Exponential Smoothing | 4749.165741 | 55.864659 | 38.516078 | 22.095183 | 179.560484 | 165.919134 |
| 5 | Deseasonal - Exponential Smoothing | 4774.382729 | 55.5275 | 37.488666 | 21.8423 | 179.560484 | 166.749778 |
| 6 | BoxCox - Exponential Smoothing | 4653.246555 | 54.946508 | 37.555881 | 21.63599 | 179.560484 | 166.093056 |
| 7 | Stationary - Exponential Smoothing | 5658.627363 | 63.86568 | 43.658278 | 26.080376 | 179.560484 | 179.444888 |
# Report do algoritmo arima, também é adicionado a base de resultados e realizado o plot de acompanhamento
def reportArima(arimaModel, modelName, calcFunction):
global result
global figs
mse = []
rmse = []
mae = []
mape = []
mtv = []
mpv = []
title = modelName + ' - Auto Arima'
indexPlot = 0
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
arimaModel.fit(train)
pred = arimaModel.predict(test.shape[0])
pred, test, indexPlot = calcFunction(train, pred, test)
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
# Correlação entre os periodos com ARIMA
lags = 90
with catch_warnings():
filterwarnings('ignore')
fig, ax = plt.subplots(2, figsize = (12, 6), sharex = True)
plot_acf(seriesOriginal.diff().dropna(), ax = ax[0], lags = lags, missing = 'drop')
plot_pacf(seriesOriginal.diff().dropna(), ax = ax[1], lags = lags)
plt.show()
%%time
# Utilizando o auto arima para descobrir os argumentos ideias baseados no conjunto de dado informado
data = seriesOriginal.copy()
arimaModel = auto_arima(seriesOriginal.copy(), m = FEATURES, seasonal = False)
arimaModel
C:\Users\herik\anaconda3\lib\site-packages\pmdarima\arima\_validation.py:62: UserWarning: m (7) set for non-seasonal fit. Setting to 0
CPU times: total: 2.28 s Wall time: 2.36 s
ARIMA(order=(1, 1, 1), scoring_args={}, suppress_warnings=True,
with_intercept=False)
reportArima(arimaModel, 'Original', calcPredTestOriginal)
%%time
data = seriesDeseasonal.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = False)
arimaModel
C:\Users\herik\anaconda3\lib\site-packages\pmdarima\arima\_validation.py:62: UserWarning: m (7) set for non-seasonal fit. Setting to 0
CPU times: total: 2.72 s Wall time: 2.78 s
ARIMA(order=(2, 1, 1), scoring_args={}, suppress_warnings=True)
reportArima(arimaModel, 'Deseasonal', calcPredTestseriesDeseasonal)
%%time
data = seriesBoxCox.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = True)
arimaModel
CPU times: total: 27.9 s Wall time: 28.2 s
ARIMA(order=(1, 1, 1), scoring_args={}, seasonal_order=(1, 0, 1, 7),
suppress_warnings=True)
reportArima(arimaModel, 'BoxCox', calcPredTestBoxCox)
%%time
data = seriesResidual.copy()
arimaModel = auto_arima(data, m = FEATURES, seasonal = True)
arimaModel
CPU times: total: 45.5 s Wall time: 46 s
ARIMA(order=(1, 0, 5), scoring_args={}, seasonal_order=(0, 0, 1, 7),
suppress_warnings=True)
reportArima(arimaModel, 'Stationary', calcPredTestStationary)
result
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | |
|---|---|---|---|---|---|---|---|
| 0 | Original - Time Series Regression | 5296.329792 | 60.784585 | 42.960631 | 25.553663 | 179.560484 | 168.008596 |
| 1 | Deseasonal - Time Series Regression | 5081.77776 | 58.208746 | 40.122247 | 23.488563 | 179.560484 | 164.670085 |
| 2 | BoxCox - Time Series Regression | 5328.602465 | 61.060147 | 43.254815 | 25.416806 | 179.560484 | 165.399796 |
| 3 | Stationary - Time Series Regression | 6057.828087 | 64.323443 | 46.538896 | 27.899678 | 179.560484 | 184.267746 |
| 4 | Original - Exponential Smoothing | 4749.165741 | 55.864659 | 38.516078 | 22.095183 | 179.560484 | 165.919134 |
| 5 | Deseasonal - Exponential Smoothing | 4774.382729 | 55.5275 | 37.488666 | 21.8423 | 179.560484 | 166.749778 |
| 6 | BoxCox - Exponential Smoothing | 4653.246555 | 54.946508 | 37.555881 | 21.63599 | 179.560484 | 166.093056 |
| 7 | Stationary - Exponential Smoothing | 5658.627363 | 63.86568 | 43.658278 | 26.080376 | 179.560484 | 179.444888 |
| 8 | Original - Auto Arima | 5707.329361 | 63.881702 | 45.981605 | 27.441516 | 179.560484 | 169.582966 |
| 9 | Deseasonal - Auto Arima | 9591.527185 | 77.052424 | 59.641758 | 35.25454 | 179.560484 | 191.567949 |
| 10 | BoxCox - Auto Arima | 11358.34545 | 84.301119 | 65.970324 | 39.351481 | 179.560484 | 206.5619 |
| 11 | Stationary - Auto Arima | 9014.243165 | 73.755288 | 56.443311 | 33.402952 | 179.560484 | 192.709473 |
# Redefinindo variaveis globais para o treino utilizando LSTM
BATCH_SIZE = 30
# hypertuning do algoritmo de LSTM
def GSLSTM(data, learning_rate, calcFunction):
mse = []
# Crossvalidation para cada parte do conjunto
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
try:
with catch_warnings():
filterwarnings('ignore')
# Normalização e reshape do conjunto de treino
ss = StandardScaler()
ss.fit(train.values.reshape(-1, 1))
train_input = ss.transform(train.values.reshape(-1, 1))
# Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
test_input = train_input[-(FEATURES + 1):]
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reset da sessão
clear_session()
# Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
initializer = GlorotUniform(seed = SEED)
model = Sequential()
# 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
model.add(Dense(1, kernel_initializer = initializer))
model.add(BatchNormalization())
model.add(Dense(1, kernel_initializer = initializer))
# Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
# Treinando o modelo com otimizador Adam
model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
pred = np.empty(test.shape[0])
# Realizando predições no conjunto de teste
for i in range(len(pred)):
prediction = model.predict(test_gen, verbose = False)
pred[i] = prediction
test_input = np.delete(test_input, 0, 0)
test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reorganizando o shape e chamando a função de calculo
pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
pred, test, _ = calcFunction(train, pred, test)
mse.append(mean_squared_error(pred, test))
except:
mse.append(-1)
return np.mean(mse)
# Report do algoritmo LSTM
def reportLSTM(data, modelName, calcFunction, learning_rate):
global result
global figs
mse = []
rmse = []
mae = []
mape = []
mtv = []
mpv = []
title = modelName + ' - Long Short Term Memory (LSTM)'
# Crossvalidation para cada parte do conjunto
for train_id, test_id in split_range.split(data):
train, test = data.iloc[train_id], data.iloc[test_id]
# Normalização e reshape do conjunto de treino
ss = StandardScaler()
ss.fit(train.values.reshape(-1, 1))
train_input = ss.transform(train.values.reshape(-1, 1))
# Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
test_input = train_input[-(FEATURES + 1):]
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reset da sessão
clear_session()
# Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
initializer = GlorotUniform(seed = SEED)
model = Sequential()
# 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
model.add(Dense(1, kernel_initializer = initializer))
model.add(BatchNormalization())
model.add(Dense(1, kernel_initializer = initializer))
# Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
# Treinando o modelo com otimizador Adam
model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
h = model.fit(train_gen, epochs = EPOCHS, callbacks = [early_stopping], verbose = False)
pred = np.empty(test.shape[0])
# Realizando predições no conjunto de teste
for i in range(len(pred)):
prediction = model.predict(test_gen, verbose = False)
pred[i] = prediction
test_input = np.delete(test_input, 0, 0)
test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reorganizando o shape e chamando a função de calculo
pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
pred, test, indexPlot = calcFunction(train, pred, test)
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_on(test, pred, mse, rmse, mae, mape, mtv, mpv)
result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mape), np.mean(mtv), np.mean(mpv), False)
return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)
%%time
if EXECUTE_GRID_SEARCH:
# Gerando tabela de hypertunning com taxas de learning_rate
learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesOriginal.copy(), x.learning_rate, calcPredTestOriginal), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 998 µs
# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()
reportLSTM(seriesOriginal.copy(), 'Original', calcPredTestOriginal, 0.0001)
%%time
if EXECUTE_GRID_SEARCH:
learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesDeseasonal.copy(), x.learning_rate, calcPredTestseriesDeseasonal), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
reportLSTM(seriesDeseasonal.copy(), 'Deseasonal', calcPredTestseriesDeseasonal, 0.01)
%%time
if EXECUTE_GRID_SEARCH:
learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesBoxCox.copy(), x.learning_rate, calcPredTestBoxCox), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 0 ns
reportLSTM(seriesBoxCox.copy(), 'BoxCox', calcPredTestBoxCox, 0.001)
%%time
if EXECUTE_GRID_SEARCH:
learning_rates = np.logspace(-5, 1, 7)
hyperparam = pd.DataFrame(learning_rates, columns = ['learning_rate'])
hyperparam['mse'] = hyperparam.swifter.apply(lambda x: GSLSTM(seriesResidual.copy(), x.learning_rate, calcPredTestStationary), axis = 1)
display(hyperparam.query('mse == mse.min() and mse != -1'))
CPU times: total: 0 ns Wall time: 1e+03 µs
reportLSTM(seriesResidual.copy(), 'Stationary', calcPredTestStationary, 0.00001)
result
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | |
|---|---|---|---|---|---|---|---|
| 0 | Original - Time Series Regression | 5296.329792 | 60.784585 | 42.960631 | 25.553663 | 179.560484 | 168.008596 |
| 1 | Deseasonal - Time Series Regression | 5081.77776 | 58.208746 | 40.122247 | 23.488563 | 179.560484 | 164.670085 |
| 2 | BoxCox - Time Series Regression | 5328.602465 | 61.060147 | 43.254815 | 25.416806 | 179.560484 | 165.399796 |
| 3 | Stationary - Time Series Regression | 6057.828087 | 64.323443 | 46.538896 | 27.899678 | 179.560484 | 184.267746 |
| 4 | Original - Exponential Smoothing | 4749.165741 | 55.864659 | 38.516078 | 22.095183 | 179.560484 | 165.919134 |
| 5 | Deseasonal - Exponential Smoothing | 4774.382729 | 55.5275 | 37.488666 | 21.8423 | 179.560484 | 166.749778 |
| 6 | BoxCox - Exponential Smoothing | 4653.246555 | 54.946508 | 37.555881 | 21.63599 | 179.560484 | 166.093056 |
| 7 | Stationary - Exponential Smoothing | 5658.627363 | 63.86568 | 43.658278 | 26.080376 | 179.560484 | 179.444888 |
| 8 | Original - Auto Arima | 5707.329361 | 63.881702 | 45.981605 | 27.441516 | 179.560484 | 169.582966 |
| 9 | Deseasonal - Auto Arima | 9591.527185 | 77.052424 | 59.641758 | 35.25454 | 179.560484 | 191.567949 |
| 10 | BoxCox - Auto Arima | 11358.34545 | 84.301119 | 65.970324 | 39.351481 | 179.560484 | 206.5619 |
| 11 | Stationary - Auto Arima | 9014.243165 | 73.755288 | 56.443311 | 33.402952 | 179.560484 | 192.709473 |
| 12 | Original - Long Short Term Memory (LSTM) | 9381.265495 | 84.147816 | 61.880612 | 38.23957 | 179.560484 | 188.430568 |
| 13 | Deseasonal - Long Short Term Memory (LSTM) | 5064.39062 | 57.410273 | 39.893818 | 23.371864 | 179.560484 | 164.653074 |
| 14 | BoxCox - Long Short Term Memory (LSTM) | 5325.842935 | 63.086576 | 44.862384 | 25.808698 | 179.560484 | 159.906651 |
| 15 | Stationary - Long Short Term Memory (LSTM) | 8713.79911 | 76.748656 | 57.541623 | 34.481319 | 179.560484 | 202.934807 |
def create_features(df, label = None):
# Separando variaveis temporais
df['date'] = df.index
df['hour'] = df['date'].dt.hour
df['dayofweek'] = df['date'].dt.dayofweek
df['quarter'] = df['date'].dt.quarter
df['month'] = df['date'].dt.month
df['year'] = df['date'].dt.year
df['dayofyear'] = df['date'].dt.dayofyear
df['dayofmonth'] = df['date'].dt.day
df['weekofyear'] = df['date'].dt.weekofyear
X = df[['hour','dayofweek','quarter','month','year',
'dayofyear','dayofmonth','weekofyear']]
if label:
y = df[label]
return X, y
return X
def reportModeloRegressao(modelo, x_teste, y_teste):
y_pred = modelo.predict(x_teste)
erros = abs(y_pred - y_teste)
mape = 100 * np.mean(erros / y_teste)
acuracia = 100 - mape
mse = mean_squared_error(y_teste, y_pred, squared = True)
mae = mean_absolute_error(y_teste, y_pred)
rmse = mean_squared_error(y_teste, y_pred, squared = False)
print(modelo,'\n')
print('Dados de teste')
print('Acuracia : {:0.2f}%'.format(acuracia))
print('MAE : {:0.2f}'.format(mae))
print('MSE : {:0.2f}'.format(mse))
print('RMSE : {:0.2f}\n'.format(rmse))
def treinaRegressao_GridSearchCV(modelo, params_, x_treino, y_treino, x_teste, y_teste,\
n_jobs = -1, cv = 5, refit = True, scoring = None, salvar_resultados = False,\
retorna_modelo = False):
grid = GridSearchCV(modelo, params_, n_jobs = n_jobs, cv = cv, refit = refit, scoring = scoring)
grid.fit(x_treino, y_treino, eval_set = [(x_treino, y_treino)], verbose = False)
pred = grid.predict(x_teste)
modelo_ = grid.best_estimator_
print(grid.best_params_)
reportModeloRegressao(modelo_, x_teste, y_teste)
if salvar_resultados:
resultados_df = pd.DataFrame(grid.cv_results_)
if retorna_modelo:
return resultados_df, modelo_
else:
resultados_df
if retorna_modelo:
return modelo_
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()
# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
%%time
if EXECUTE_GRID_SEARCH:
'''params = {
'max_depth': np.arange(4, 8, step = 1),
'gamma': np.arange(1.1, step = 0.1),
'reg_alpha' : np.arange(1.1, step = 0.1),
'reg_lambda' : np.arange(1.1, step = 0.1),
'colsample_bytree' : np.arange(1.1, step = 0.1),
'min_child_weight' : np.arange(11, step = 1),
'eta': np.arange(1.1, step = 0.1),
'n_estimators': [1000],
'seed': [SEED],
'early_stopping_rounds': [10],
'objective': ['reg:squarederror']
}'''
params = {
'max_depth': np.arange(4, 8, step = 1),
'gamma': np.arange(1.1, step = 0.1),
'reg_alpha' : np.arange(1.1, step = 0.1),
'reg_lambda' : np.arange(1.1, step = 0.1),
'eta': np.arange(1.1, step = 0.1),
'n_estimators': [1000],
'seed': [SEED],
'early_stopping_rounds': [10],
'objective': ['reg:squarederror']
}
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', seed = SEED)
treinaRegressao_GridSearchCV(model, params, X_train, y_train, X_test, y_test, scoring = 'neg_root_mean_squared_error')
CPU times: total: 0 ns Wall time: 2 ms
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
eval_set = [(X_train, y_train)],
verbose = False)
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
early_stopping_rounds=10, enable_categorical=False, eta=0.1,
eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
importance_type=None, interaction_constraints='',
learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
missing=nan, monotone_constraints='()', n_estimators=1000,
n_jobs=0, num_parallel_tree=1, predictor='auto',
random_state=84796315, reg_alpha=0.7, ...)
plot_importance(model, height = 1)
<AxesSubplot:title={'center':'Feature importance'}, xlabel='F score', ylabel='Features'>
test['Previsões'] = model.predict(X_test)
XGBresults = pd.concat([test, train], sort = False)
pred = test['Previsões']
test_ = test['Vendas']
mse, rmse, mae, mape, mtv, mpv = calculate_metrics_off(test_, pred)
result = record(result, 'Original - XGBoost', mse, rmse, mae, mape, mtv, mpv, False)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
y = seriesOriginal,
name = 'Real'),
go.Scatter(x = test.index,
y = test['Previsões'],
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Original - XGBoost')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
pyoff.iplot(fig)
result
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | |
|---|---|---|---|---|---|---|---|
| 0 | Original - Time Series Regression | 5296.329792 | 60.784585 | 42.960631 | 25.553663 | 179.560484 | 168.008596 |
| 1 | Deseasonal - Time Series Regression | 5081.77776 | 58.208746 | 40.122247 | 23.488563 | 179.560484 | 164.670085 |
| 2 | BoxCox - Time Series Regression | 5328.602465 | 61.060147 | 43.254815 | 25.416806 | 179.560484 | 165.399796 |
| 3 | Stationary - Time Series Regression | 6057.828087 | 64.323443 | 46.538896 | 27.899678 | 179.560484 | 184.267746 |
| 4 | Original - Exponential Smoothing | 4749.165741 | 55.864659 | 38.516078 | 22.095183 | 179.560484 | 165.919134 |
| 5 | Deseasonal - Exponential Smoothing | 4774.382729 | 55.5275 | 37.488666 | 21.8423 | 179.560484 | 166.749778 |
| 6 | BoxCox - Exponential Smoothing | 4653.246555 | 54.946508 | 37.555881 | 21.63599 | 179.560484 | 166.093056 |
| 7 | Stationary - Exponential Smoothing | 5658.627363 | 63.86568 | 43.658278 | 26.080376 | 179.560484 | 179.444888 |
| 8 | Original - Auto Arima | 5707.329361 | 63.881702 | 45.981605 | 27.441516 | 179.560484 | 169.582966 |
| 9 | Deseasonal - Auto Arima | 9591.527185 | 77.052424 | 59.641758 | 35.25454 | 179.560484 | 191.567949 |
| 10 | BoxCox - Auto Arima | 11358.34545 | 84.301119 | 65.970324 | 39.351481 | 179.560484 | 206.5619 |
| 11 | Stationary - Auto Arima | 9014.243165 | 73.755288 | 56.443311 | 33.402952 | 179.560484 | 192.709473 |
| 12 | Original - Long Short Term Memory (LSTM) | 9381.265495 | 84.147816 | 61.880612 | 38.23957 | 179.560484 | 188.430568 |
| 13 | Deseasonal - Long Short Term Memory (LSTM) | 5064.39062 | 57.410273 | 39.893818 | 23.371864 | 179.560484 | 164.653074 |
| 14 | BoxCox - Long Short Term Memory (LSTM) | 5325.842935 | 63.086576 | 44.862384 | 25.808698 | 179.560484 | 159.906651 |
| 15 | Stationary - Long Short Term Memory (LSTM) | 8713.79911 | 76.748656 | 57.541623 | 34.481319 | 179.560484 | 202.934807 |
| 16 | Original - XGBoost | 4667.546621 | 68.319445 | 59.339327 | 28.327163 | 221.854839 | 190.8423 |
result.to_csv('export\resultados_time_series.csv')
# Tratando nomes e criando colunas de controle para os resultados gerados
topResult = (
result
.assign(Full_Name = lambda x: x.Algorithm.apply(lambda a: a.split('(')[0]
.replace('ARIMA', 'Auto Arima')
.replace('Long Short Term Memory', 'LSTM')))
.assign(Data_Category = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[0]))
.assign(Algorithm = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[1].split('(')[0]
.replace('ARIMA', 'Auto Arima')
.replace('Long Short Term Memory', 'LSTM')))
.sort_values('MSE')
)
topResult.head()
| Algorithm | MSE | RMSE | MAE | MAPE | Mean_True_Value | Mean_Predict_Value | Full_Name | Data_Category | |
|---|---|---|---|---|---|---|---|---|---|
| 6 | Exponential Smoothing | 4653.246555 | 54.946508 | 37.555881 | 21.63599 | 179.560484 | 166.093056 | BoxCox - Exponential Smoothing | BoxCox |
| 16 | XGBoost | 4667.546621 | 68.319445 | 59.339327 | 28.327163 | 221.854839 | 190.8423 | Original - XGBoost | Original |
| 4 | Exponential Smoothing | 4749.165741 | 55.864659 | 38.516078 | 22.095183 | 179.560484 | 165.919134 | Original - Exponential Smoothing | Original |
| 5 | Exponential Smoothing | 4774.382729 | 55.5275 | 37.488666 | 21.8423 | 179.560484 | 166.749778 | Deseasonal - Exponential Smoothing | Deseasonal |
| 13 | LSTM | 5064.39062 | 57.410273 | 39.893818 | 23.371864 | 179.560484 | 164.653074 | Deseasonal - LSTM | Deseasonal |
# Plot dos resultados obtidos por ordem ascendente do MSE
colors = {'Time Series Regression':'#1688F2',
'Exponential Smoothing':'#17FC2E',
'Auto Arima': '#E6C220',
'LSTM ': '#FC4417',
'XGBoost': '#A13BF5'}
# plotly figure
fig = go.Figure(layout = go.Layout(yaxis = {'type': 'category', 'title': 'Algoritmo e Categoria'},
xaxis = {'title': 'MSE'},
title = 'MSE por Algoritmo e Tipo de Dado'))
for t in topResult['Algorithm'].unique():
topResultFiltered = topResult[topResult['Algorithm']== t].copy()
fig.add_traces(go.Bar(x = topResultFiltered['MSE'], y = topResultFiltered['Full_Name'], name = str(t),\
marker_color = str(colors[t]), orientation = 'h',
text = round(topResultFiltered['MSE'].astype(np.double)), textposition = 'outside'))
fig.update_layout(yaxis = {'categoryorder':'total descending'}, autosize = False,
width = 1450,
height = 800)
fig.show()
pred_range = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 10, 31))
split_range = TimeSeriesSplit(n_splits = 2, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])
BATCH_SIZE = 1000
# Alocando melhor modelo a memória e separando base de treino
data = seriesBoxCox.copy()
train = data[datetime.date(2017, 1, 1): datetime.date(2018, 6, 16)]
# Treinando modelo baseado dos parâmetros descobertos na fase de modelagem
ES = (
ExponentialSmoothing(train, trend = 'add', seasonal = 'add', seasonal_periods = FEATURES, damped_trend = True)
.fit(smoothing_level = 0.1, smoothing_trend = 0.7, smoothing_seasonal = 0.0, method = 'ls', damping_trend = 0.8)
)
# Calculando a previsão até o final do ano de 2018
pred = ES.predict(str(data.index[0]), '2018-12-31')
pred, data, _ = calcPredTestBoxCox(train, pred, data)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
y = seriesOriginal,
name = 'Real'),
go.Scatter(x = pred.index,
y = pred,
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'BoxCox - Exponential Smoothing')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31',
annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'green', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
ES.summary()
| Dep. Variable: | None | No. Observations: | 532 |
|---|---|---|---|
| Model: | ExponentialSmoothing | SSE | 7800.497 |
| Optimized: | True | AIC | 1452.579 |
| Trend: | Additive | BIC | 1503.899 |
| Seasonal: | Additive | AICC | 1453.392 |
| Seasonal Periods: | 7 | Date: | Sun, 23 Oct 2022 |
| Box-Cox: | False | Time: | 13:02:16 |
| Box-Cox Coeff.: | None |
| coeff | code | optimized | |
|---|---|---|---|
| smoothing_level | 0.1000000 | alpha | False |
| smoothing_trend | 0.7000000 | beta | False |
| smoothing_seasonal | 0.000000 | gamma | False |
| initial_level | 3.5625210 | l.0 | True |
| initial_trend | 1.3021796 | b.0 | True |
| damping_trend | 0.8000000 | phi | False |
| initial_seasons.0 | -8.0164149 | s.0 | True |
| initial_seasons.1 | -3.7667775 | s.1 | True |
| initial_seasons.2 | -3.8599627 | s.2 | True |
| initial_seasons.3 | -4.3632835 | s.3 | True |
| initial_seasons.4 | -5.0029325 | s.4 | True |
| initial_seasons.5 | -6.0412289 | s.5 | True |
| initial_seasons.6 | -9.3733856 | s.6 | True |
BATCH_SIZE = 30
# Definindo dataset
dtHistoryLSTM = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date
# Separando periodo minimo e maximo
startLSTM = dtHistoryLSTM.min()
endLSTM = dtHistoryLSTM.max()
# Separando IDs das lacunas
idxLSTM = pd.date_range(startLSTM, endLSTM, normalize = True)
# Transformando dt para série e realizando contagem de valores diários
seriesOriginalLSTM = dtHistoryLSTM.value_counts(sort = False).sort_index().reindex(idxLSTM, fill_value = 0)
# Removendo outliers
seriesOriginalLSTM = seriesOriginalLSTM[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
# Adicionando predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalTSTRPredictions = seriesOriginalLSTM.append(newPredictions)
# Gerando STL para separar o Deseasonal
stl = STL(seriesOriginalTSTRPredictions)
res = stl.fit()
# Separando Deseasonal
seriesOriginalTSTRPredictionsDeseasonal = res.observed - res.seasonal
# Separando index da série
dataTime = seriesOriginalTSTRPredictionsDeseasonal.copy().index
# Separando base de treino e teste
train_size = 532
test_size = len(seriesOriginalTSTRPredictionsDeseasonal) - train_size
train = seriesOriginalTSTRPredictionsDeseasonal[0:train_size]
test = seriesOriginalTSTRPredictionsDeseasonal[train_size:len(seriesOriginalTSTRPredictionsDeseasonal)]
# Normalização e reshape do conjunto de treino
scaler = MinMaxScaler(feature_range=(0, 1))
scaler.fit(train.values.reshape(-1, 1))
train_input = scaler.transform(train.values.reshape(-1, 1))
# Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
test_input = train_input[-(FEATURES + 1):]
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
train_gen = TimeseriesGenerator(train_input, train_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reset da sessão
clear_session()
# Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
initializer = GlorotUniform(seed = SEED)
model = Sequential()
# 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
model.add(LSTM(64, input_shape = (FEATURES, 1), return_sequences = False))
model.add(Dense(1, kernel_initializer = initializer))
model.add(BatchNormalization())
model.add(Dense(1, kernel_initializer = initializer))
# Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
# Treinando o modelo com otimizador Adam
model.compile(loss = 'mse', optimizer = Adam(learning_rate = 0.01), metrics = ['mae'])
model.fit(train_gen, epochs = 100, callbacks = [early_stopping], verbose = True)
pred = np.empty(test.shape[0])
# Realizando predições no conjunto de teste
for i in range(len(pred)):
prediction = model.predict(test_gen, verbose = False)
pred[i] = prediction
#print(test_gen[0][0], ' = ', prediction)
test_input = np.delete(test_input, 0, 0)
test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
test_gen = TimeseriesGenerator(test_input, test_input, length = FEATURES, batch_size = BATCH_SIZE)
# Reorganizando o shape e chamando a função de calculo
pred = scaler.inverse_transform(pred.reshape(-1,1)).reshape(-1)
pred, test, _ = calcPredTestseriesDeseasonal(train, pred, test)
Epoch 1/100 18/18 [==============================] - 2s 5ms/step - loss: 0.0131 - mae: 0.0827 Epoch 2/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0096 - mae: 0.0789 Epoch 3/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0071 - mae: 0.0633 Epoch 4/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0069 - mae: 0.0667 Epoch 5/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0067 - mae: 0.0657 Epoch 6/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0067 - mae: 0.0655 Epoch 7/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0065 - mae: 0.0645 Epoch 8/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0072 - mae: 0.0656 Epoch 9/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0069 - mae: 0.0670 Epoch 10/100 18/18 [==============================] - 0s 4ms/step - loss: 0.0073 - mae: 0.0670
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = dataTime,
y = seriesOriginalTSTRPredictionsDeseasonal[:594],
name = 'Real'),
go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
y = pred,
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Deseasonal - LSTM')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31',
annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'green', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
BATCH_SIZE = 1000
# Definindo dataset
dtHistoryTSR = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date
# Separando periodo minimo e maximo
startTSR = dtHistoryTSR.min()
endTSR = dtHistoryTSR.max()
# Separando IDs das lacunas
idxTSR = pd.date_range(startTSR, endTSR, normalize = True)
# Transformando dt para série e realizando contagem de valores diários
seriesOriginalTSR = dtHistoryTSR.value_counts(sort = False).sort_index().reindex(idxTSR, fill_value = 0)
# Removendo outliers
seriesOriginalTSR = seriesOriginalTSR[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
# Adicionando predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalTSRRPredictions = seriesOriginalTSR.append(newPredictions)
# Gerando STL para separar o Deseasonal
stl = STL(seriesOriginalTSRRPredictions)
res = stl.fit()
# Separando Deseasonal
seriesOriginalTSRRPredictionsDeseasonal = res.observed - res.seasonal
# Separando base de treino e teste
train_size = 532
test_size = len(seriesOriginalTSRRPredictionsDeseasonal) - train_size
train = seriesOriginalTSRRPredictionsDeseasonal[0:train_size]
test = seriesOriginalTSRRPredictionsDeseasonal[train_size:len(seriesOriginalTSRRPredictionsDeseasonal)]
# Gerando generator de treino
gen = TimeseriesGenerator(train, train, FEATURES, batch_size = BATCH_SIZE)
X_train = gen[0][0]
y_train = gen[0][1]
# Treinando modelo
lr = LinearRegression()
lr.fit(X_train, y_train)
# Separando primeira leva de treino
X_pred = y_train[-FEATURES:].reshape(1,-1)
pred = np.empty(test.shape[0])
# Realizando predições e realocando vetor de entrada
for i in range(len(pred)):
forecast = lr.predict(X_pred)
X_pred = np.delete(X_pred, 0, 1)
X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)
pred[i] = forecast
# Convertendo deseasonal de volta para os valores originais
pred, test, _ = calcPredTestseriesDeseasonal(train, pred, test)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = dataTime,
y = seriesOriginalTSTRPredictionsDeseasonal[:594],
name = 'Real'),
go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
y = pred,
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Deseasonal - Time Series Regressor')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31',
annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'green', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
# Definindo dataset
dtHistoryArima = pd.to_datetime(dtOrdersAdjusted['order_purchase_timestamp']).dt.date
# Separando periodo minimo e maximo
startArima = dtHistoryArima.min()
endArima = dtHistoryArima.max()
# Separando IDs das lacunas
idxArima = pd.date_range(startArima, endArima, normalize = True)
# Transformando dt para série e realizando contagem de valores diários
seriesOriginalArima = dtHistoryArima.value_counts(sort = False).sort_index().reindex(idxArima, fill_value = 0)
# Removendo outliers
seriesOriginalArima = seriesOriginalArima[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
# Adicionando predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
seriesOriginalArimaPredictions = seriesOriginalArima.append(newPredictions)
# Gerando STL para separar cada um dos componentes
stl = STL(seriesOriginalArimaPredictions)
res = stl.fit()
# Separando boxcox
seriesBoxCoxArima, lmbda = boxcox(seriesOriginalArimaPredictions + 1)
seriesBoxCoxArima = pd.Series(seriesBoxCoxArima, index = seriesOriginalArimaPredictions.index)
# Separando base de treino e teste
train_size = 532
test_size = len(seriesBoxCoxArima) - train_size
train = seriesBoxCoxArima[0:train_size]
test = seriesBoxCoxArima[train_size:len(seriesBoxCoxArima)]
arimaModel = auto_arima(seriesBoxCoxArima, m = FEATURES, seasonal = True)
arimaModel
ARIMA(order=(1, 1, 1), scoring_args={}, seasonal_order=(1, 0, 2, 7),
suppress_warnings=True, with_intercept=False)
arimaModel.fit(train)
pred = arimaModel.predict(test.shape[0])
pred, test, _ = calcPredTestBoxCox(train, pred, test)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = pd.date_range(datetime.date(2017, 1, 1), datetime.date(2018, 8, 17)),
y = seriesOriginalArima[:594],
name = 'Real'),
go.Scatter(x = pd.date_range(datetime.date(2018, 6, 17), datetime.date(2018, 12, 31)),
y = pred,
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'BoxCox - Arima')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31',
annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'green', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()
# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
eval_set = [(X_train, y_train)],
verbose = False)
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
early_stopping_rounds=10, enable_categorical=False, eta=0.1,
eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
importance_type=None, interaction_constraints='',
learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
missing=nan, monotone_constraints='()', n_estimators=1000,
n_jobs=0, num_parallel_tree=1, predictor='auto',
random_state=84796315, reg_alpha=0.7, ...)
# Adicionando predições futuras
newPredictions = pd.Series([0 for i in range(136)])
newPredictions.index = pd.date_range(datetime.date(2018, 8, 18), datetime.date(2018, 12, 31))
newPredictions = pd.DataFrame(newPredictions)
newPredictions = create_features(newPredictions, label = None)
newPredictionsTest = pd.concat([X_test, newPredictions])
test = pd.concat([test, newPredictions])
test['Previsões'] = model.predict(newPredictionsTest)
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = seriesOriginal.index,
y = seriesOriginal,
name = 'Real'),
go.Scatter(x = test.index,
y = test['Previsões'],
name = 'Previsto')]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Original - XGBoost')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.add_vrect(x0 = '2018-08-17', x1 = '2018-12-31',
annotation_text = 'Projeção de<br>Vendas Futuras', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'green', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
# Convertendo dtHistory para serial
dtHistoryIndexado = dtHistory.copy()
dtHistoryIndexado.index = dtHistoryIndexado['Data']
dtHistoryIndexado = dtHistoryIndexado['Vendas']
# Separando base de treino e teste
split_date = '2018-06-16'
train = dtHistoryIndexado.loc[dtHistoryIndexado.index <= split_date].copy()
test = dtHistoryIndexado.loc[dtHistoryIndexado.index > split_date].copy()
# Convertendo base de treio e teste para dataframe
train = pd.DataFrame(train)
test = pd.DataFrame(test)
# Separando X e y em treino e teste
X_train, y_train = create_features(train, label = 'Vendas')
X_test, y_test = create_features(test, label = 'Vendas')
# Treinando modelo
model = XGBRegressor(n_estimators = 1000, early_stopping_rounds = 10, objective = 'reg:squarederror', max_depth = 4, gamma = 0.7, reg_alpha = 0.7, reg_lambda = 0.9, eta = 0.1, seed = SEED)
model.fit(X_train, y_train,
eval_set = [(X_train, y_train)],
verbose = False)
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
early_stopping_rounds=10, enable_categorical=False, eta=0.1,
eval_metric=None, gamma=0.7, gpu_id=-1, grow_policy='depthwise',
importance_type=None, interaction_constraints='',
learning_rate=0.100000001, max_bin=256, max_cat_to_onehot=4,
max_delta_step=0, max_depth=4, max_leaves=0, min_child_weight=1,
missing=nan, monotone_constraints='()', n_estimators=1000,
n_jobs=0, num_parallel_tree=1, predictor='auto',
random_state=84796315, reg_alpha=0.7, ...)
newPredictionsTest = pd.concat([X_train, X_test])
test = pd.concat([train, test])
test['Previsões'] = model.predict(newPredictionsTest)
test['ci'] = 0.04 * test['Vendas']
test['ci_lower'] = test['Vendas'] - test['ci']
test['ci_upper'] = test['Vendas'] + test['ci']
test.tail()
| Vendas | date | hour | dayofweek | quarter | month | year | dayofyear | dayofmonth | weekofyear | Previsões | ci | ci_lower | ci_upper | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Data | ||||||||||||||
| 2018-08-13 | 284 | 2018-08-13 | 0 | 0 | 3 | 8 | 2018 | 225 | 13 | 33 | 239.025238 | 11.36 | 272.64 | 295.36 |
| 2018-08-14 | 311 | 2018-08-14 | 0 | 1 | 3 | 8 | 2018 | 226 | 14 | 33 | 231.174133 | 12.44 | 298.56 | 323.44 |
| 2018-08-15 | 282 | 2018-08-15 | 0 | 2 | 3 | 8 | 2018 | 227 | 15 | 33 | 255.116379 | 11.28 | 270.72 | 293.28 |
| 2018-08-16 | 316 | 2018-08-16 | 0 | 3 | 3 | 8 | 2018 | 228 | 16 | 33 | 237.553925 | 12.64 | 303.36 | 328.64 |
| 2018-08-17 | 249 | 2018-08-17 | 0 | 4 | 3 | 8 | 2018 | 229 | 17 | 33 | 203.458313 | 9.96 | 239.04 | 258.96 |
# Plot
# Definição dos dados no plot
plot_data = [go.Scatter(x = test.index,
y = test['Previsões'],
name = 'Previsto'),
go.Scatter(
name = '95% CI Upper',
x = test.index,
y = test['ci_upper'],
mode = 'lines',
marker = dict(color = '#444'),
line = dict(width = 0),
showlegend = False
),
go.Scatter(
name = '95% CI Lower',
x = test.index,
y = test['ci_upper'],
marker = dict(color = '#444'),
line = dict(width = 0),
mode = 'lines',
fillcolor = 'rgba(68, 68, 68, 0.3)',
fill = 'tonexty',
showlegend = False
)
]
# Layout
plot_layout = go.Layout(xaxis = {'title': 'Período'},
yaxis = {'title': 'Vendas'},
title = 'Original - XGBoost')
# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
pyoff.iplot(fig)
fig = go.Figure([
go.Scatter(
name = 'Previsto',
x = test.index,
y = test['Previsões'],
),
go.Scatter(
name = '95% CI Upper',
x = test.index,
y = test['ci_upper'],
mode = 'lines',
marker = dict(color = '#444'),
line = dict(width = 0),
showlegend = True
),
go.Scatter(
name = '95% CI Lower',
x = test.index,
y = test['ci_lower'],
marker = dict(color = '#444'),
line = dict(width = 0),
mode = 'lines',
fillcolor = 'rgba(255, 0, 0, 0.80)',
fill = 'tonexty',
showlegend = True
)
])
fig.update_layout(
xaxis_title = 'Período',
yaxis_title = 'Vendas',
title = 'Original - XGBoost',
hovermode = 'x',
yaxis_range = [0 , 800]
)
fig.add_vrect(x0 = '2018-06-17', x1 = '2018-08-17',
annotation_text = 'Previsão base<br>de teste', annotation_position = 'top left',
annotation = dict(font_size = 23, font_family = 'Times New Roman'),
fillcolor = 'orange', opacity = 0.2, line_width = 0)
fig.update_layout(yaxis_range = [0 , 800])
fig.update_yaxes(rangemode='tozero')